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ABSTRACT 

In compressed sensing one measures sparse signals directly 
in a compressed form via a linear transform and then recon- 
structs the original signal. However, it is often the case that 
the linear transform itself is known only approximately, a sit- 
uation called matrix uncertainty, and that the measurement 
process is noisy. Here we present two contributions to this 
problem: first, we use the replica method to determine the 
mean-squared error of the Bayes-optimal reconstruction of 
sparse signals under matrix uncertainty. Second, we consider 
a robust variant of the approximate message passing algo- 
rithm and demonstrate numerically that in the limit of large 
systems, this algorithm matches the optimal performance in a 
large region of parameters. 

Index Terms — Compressed sensing. Measurement un- 
certainty. Belief propagation. Message passing. Performance 
analysis. 

1. INTRODUCTION 

Compressed sensing (CS) is designed to measure sparse sig- 
nals directly in a compressed form by acquiring a small (with 
respect to the dimension of the signal) number of random lin- 
ear projections of the signal, and subsequently reconstructing 
computationally the signal. In particular, it has been shown 
im 111 that this reconstruction is possible and computation- 
ally feasible in many cases using £i -minimization based al- 
gorithms. It is often the case in practical applications that 
the linear transform itself is known only approximately; for 
instance this is the case if calibration has not been done per- 
fectly. Several works have analyzed the performance of £i- 
regularization based algorithm under matrix uncertainty, e.g. 
f3 4 1 and references therein. In this paper we shall use instead 
a Bayesian approach. In this framework CS under matrix un- 
certainty was studied recently by |5 1. In particular the authors 
of Is) propose a generalization of the approximate message 



* The research leading to these results has received funding from the Eu- 
ropean Research Council under the European Union's 7*'' Framework Pro- 
gramme (FP/2007-2013)/ERC Grant Agreement 307087-SPARCS. 



passing algorithm (AMP) JS] H [U that treats matrix uncer- 
tainty (MU-AMP) and outperforms other methods. Empiri- 
cally, they show that the MU-AMP algorithm performs near 
oracle bounds. 

Our goal in this paper is twofold. First we compute 
the best theoretically possible performance of reconstruc- 
tion algorithms under matrix uncertainty using the replica 
method |j9l|T0l[n][l2l|8l[T3l|, and show that in a large region 
of parameters this performance is indeed matched by MU- 
AMP. And second we consider a slightly different algorithm 
that we refer to as robust- AMP algorithm (using a minimal 
change first suggested in [61) and show that it is asymptoti- 
cally equivalent to MU-AMP and thus also it reaches Bayes 
optimal MSE while not assuming the knowledge of the mea- 
surement noise, nor the level of matrix uncertainty. 

2. DEFINITIONS 

In our analysis we consider a sparse TV-dimensional signal x 
having on average K non-zero iid components chosen from a 
distribution 4){x). Defining the density p ~ K/N we have 
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Note that our analysis applies as well to non-iid signals with 
empirical distribution of components converging to (t){x), as 
discussed e.g. in [14|. For concreteness, here we shall use a 
Gaussian 4){x) of zero mean and unit variance. Although our 
results depend quantitatively on the form of (i>{x), the over- 
all qualitative picture is robust with respect to this choice. 
We further assume that the parameters of P(x) are known 
(if this is not the case they can be learned via expectation- 
maximization Is] [T3] [TSll ). One could use an approximately 
sparse signal as in lfT6l as well. 

The M measurements are obtained via noisy linear 
random projections of the signal 
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where we model the noise as a Gaussian random variable 
with mean and variance A, and where F° has iid random 
components of mean and variance 1/iV. The parameter a — 
M/N is the measurement, or sampling, rate. 

Additionally, we consider that the matrix F" is not known 
perfectly, and that we know only its noisy version; 
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where is a white noise with mean and variance 
and where 77 measures the uncertainty on the matrix. When 
77 0, the measurement matrix is perfectly known (this is 
the usual compressed sensing situation) while for i] 00 
nothing is known about the measurement matrix. Assuming 
a Gaussian distribution for the matrix elements and the noise 
X, we obtain the posterior probability of one matrix element 
P{F^i\F^^) to be a Gaussian with mean and variance 
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Our goal is to reconstruct the signal s based on the knowl- 
edge of y and F'. The Bayes-optimal estimate of a signal 
X*, that minimizes the mean-squared error (MSE) MSE = 



/N with the original signal s, is given by 



dxiXiVi{xi), (5) 



where i^i {xi ) is the marginal distribution of Xi with respect to: 
P(x|y,F') = ^ J dF°P(F°|F')e-5^(''-^'''')' . (6) 

3. REPLICA METHOD: OPTIMUM 
RECONSTRUCTION BOUNDS 

The MSE obtained from the Bayes-optimal approach can be 
computed in the limit of large N via the replica method. Al- 
though this method is in general non-rigorous, it is sometimes 
possibles to prove that the results derived from it are exact, 
see 1 9] [H [mini mill. We shall leave out a detailed deriva- 
tion and refer instead to |l8][T3] for a very similar computation. 
The result is that in the large signal limit, the MSE is given by 



MSE(a, p, A, 77) = argmax $(^) 

E 

where the "potential" ^{E) is given by 
log[A + E+{p- E)D]-i 
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Fig. 1. The potential ^{E), eq. (HJ, for measurement rate 
a ~ 0.5, matrix uncertainty r] — 10^"', and noise A — 10^^", 
for different values of density p. The global maximum $(£'*) 
determines the values of the Bayes-optimal MSE = E*. As 
the density varies the value E* (p) has a discontinuity which is 
usually referred to as a "first-order" or "discontinuous" phase 
transition. AMP algorithms are performing a steepest ascent 
of the potential and get trapped at the maximum with the 
largest value of E, instead of reaching the global one (see e.g. 
the case p = 0.33). The appearance of this trapping local 
maxima is referred to as the "spinodal" transition. 
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Note how closely related is this expression to the one derived 
in lfT3l for the usual compressed sensing with known mea- 
surement matrix (which is recovered taking 77 = 0). 

The potential ^ is shown in Fig.[T] Maximizing the po- 
tential, one obtains the Bayes-optimal MSE which we show 
in Fig.|2]as a function of 77 for different densities p. We see a 
phenomenology similar to that described in |.13i.l6J for other 
types of noise. Because of a two-maxima shape of the func- 
tion ^{E), the Bayes-optimal MSE displays a sharp transition 
separating a region of parameters with a small MSE, compa- 
rable to max(A,77), from a region with a large 0(1) value 
of the MSE. In this second region of parameters compressed 
sensing technics -regardless of the reconstruction method- 
will not be useful. This is a "first order" phase transition. 

Another important transition (the "spinodal") that can be 
studied from the form of the potential function ^{E) is the 
appearance of a high-MSE local maxima. Since the AMP al- 
gorithm performs a steepest ascent in the ^{E) function (see 
e.g. 1 13 1) this transition separates the region of parameters in 
which the performance of AMP will be asymptotically opti- 
mal from the one where AMP is suboptimal. An important 
feature of this transition, shown in Fig. |3] is that its position 
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Fig. 2. Bayes-optimal MSE for a = 0.5, A = 10"^" as a 
function of matrix uncertainty -q for different densities p. The 
MSE of the trapping local maximum is shown in dashed line. 
The results of the AMP algorithm on some instances are also 
shown (points, done with N = 20000). The results of these 
AMP runs (points) agree perfectly with the prediction (lines). 
In the inset, we show the location of the two phase transitions 
in the p, rj plane: the upper Une determines the "first-order" 
threshold beyond which the optimal-MSE suddenly degrades, 
while the lower line is the "spinodal" transition that marks the 
degradation of the performance of AMP. Note that the spin- 
odal is almost independent of the noise 77 in a large range of 
values. 



is only very weakly dependent on the value of the matrix un- 
certainty. This means that the performance of AMP is ex- 
tremely robust to matrix uncertainty 77. The same result was 
also obtained for additive measurement noise in [T3l and for 
approximate sparsity in l,16J . 



4. ROBUST MESSAGE-PASSING ALGORITHM 

The Bayesian approach to compressed sensing combined with 
belief propagation based reconstruction algorithm leads to the 
so-called approximate message passing (AMP) algorithm as 
first derived in [6] for the minimization of li, and subse- 
quently generalized in IfTSl lTl fm . This approach was adapted 
by Parker, Cevher and Schniter to treat the matrix uncertainty 
MU-AMP in 15). Here we shall consider a version of the 
canonical AMP, that we call robust-AMP, that turns out to 
be robust to noisy measurement and matrix uncertainty, so 
that it can be used indifferently of the presence or absence 
of noise and matrix uncertainty. We show in the next section 
that robust-AMP and MU-AMP are equivalent in the limit of 
infinite systems. 

For every measurement component we define one real 
number o;^, for each signal component i we define four real 
numbers Y^i, Ri, Ui, Vi. These quantities are updated as fol- 



Fig. 3. The phase diagram "a la" Donoho-Tanner ifTTl in the 
limit of large signals, — > 00. The horizontal axis is the 
sampling ratio M/N and the vertical one measures the gain 
K/M provided by compressed sensing. The dotted black line 
is the transition for the £1 reconstruction IITtI in zero noise. 
The blue lines (from |8 13 1) show the location of the spin- 
odal (dotted blue) and first-order (full blue) phase transitions 
for the noiseless Bayesian approach. While a perfect recon- 
struction (MSE=0) can be obtained in principle until p = a, 
the AMP algorithm allows such perfect reconstruction only 
up to the spinodal line (dotted blue). The noisy counterpart 
of these two transitions is shown in red, using A = 10^^ 
and 77 = 10^^. While the region where a good reconstruc- 
tion is possible (below the full red line) has shrunk a lot, the 
spinodal line (dotted red line) that marks the onset of good 
performance of AMP is left almost unchanged below p w 0.8 
demonstrating the robustness of the AMP to noises. 



lows (for the derivation with these notations see |fT3l ): 
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Where only the functions fa and depend explicitly on the 
model P{x): 

/a(S2,i?) = j AxxM{TP-,R,x) , (16) 
fc{'^\R) = J dxx'MiY.^R^x)- fl{Y.\R), (17) 



where M (S^, i?, x) is the following probability distribution 
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The explicit expression for /„ and fc for the Gauss-Bernoulli 
signal is given in llT3i while the case of approximate sparsity 
was considered in ||T6l . The above equations are initialized 
with a*^*^ = 0, w*^" — p, ioj^'^ = y^, then the equations are 
iterated till convergence. 

The difference between the present algorithm with respect 
to the more common version of AMP is in the way the esti- 
mate of the current error on the measurement element is com- 
puted in eq. (fTOl i. Most previous works (e.g. Q [T3j) used 
a /i-dependent vector = A + F^^l in the case of 

noisy measurement with a perfectly known matrix, whereas 
the original paper [6| used a value precomputed by the state 
evolution. The modification done by the authors of ref. JS] 
to incorporate the effect of matrix uncertainty is (in our nota- 
tion) = A + Fl^vl + + aj)ij/[N{l + 77)]. As 

shown in ref. jS | this leads to a very efficient algorithm with 
matrix uncertainty. The expression we use instead in eq. (fTOl i 
was also proposed in |6|. Perhaps surprisingly, it is equiva- 
lent to the one of ref. |5 1 when the system size N ^ 00 even 
in the case of matrix uncertainty. The advantage of expres- 
sion ( fTO)i is that it is not need to explicitly know the value of 
77 and the algorithm automatically incorporates the errors 
coming from the uncertainty on the matrix and the measure- 
ment; and hence we are using the same code regardless the 
presence or absence of noise and/or matrix uncertainty. We 
thus refer to algorithm ([T0)-(1T5) as the "robust"- AMP. 

5. DENSITY EVOLUTION 

The AMP approach is amenable to asymptotic {N — > 00) 
analysis in the case of i.i.d. random measurement matrices 
using a method known as the "cavity method" (in statistical 
physics) [91, the "density evolution" in coding [19], and the 
"state evolution" in the context of CS |l6l|20l- We shall not 
detail the computation here, it goes in the same lines as in 
lITSi . Given the parameters p, a, 77, A, the MSE follows: 



(19) 



where tti* follows eq. with on the right hand side, 
i?*"^" = p, Vz is a Gaussian integral, and fc is defined by 
eq. iTl\ . The comparison between the evolution of the MSE 
in the algorithm and eq. iT% is shown on Fig. |4l the agree- 
ment is very good. Note that we have also applied the density 
evolution to the AMP algorithm with matrix uncertainty of 
f5|, and found that it obeys asymptotically the same equation 
(fT9] l. More specifically, the term V* from eq. (fTOl i evolves as 
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' Of course A and rj can be learned with expectation maximization within 
the MU-AMP [8l[l3l[l5), but this adds considerable computational time. 



Fig. 4. Comparison of the time evolution of MSE computed 
with density evolution and the one found numerically using 
robust-AMP for systems of size N = 25000, p ^ 0.1. The 
agreement between theoretical predictions (eq. (fT9t . full line) 
with the data from the AMP algorithm (points) is algorithm. 
The arrows show the MSE obtained using £1 minimization on 
the same instances, for comparison. The better performance 
of the AMP approach over ii is clear. 



A + i?* + (p - E*^)D in the hmit iV ^ 00 both for the ex- 
pression in ref. [5| and in eq. (fTOl i. From eq. (fT9t , one can 
also derive that the evolution is equivalent to a steepest as- 
cent of the potential (f){E) obtained from the replica method 
in eq. This underlines the importance of the spinodal tran- 
sition illustrated on Figs.[T]|2] and|3] In particular we see that 
the region where the AMP converges to the Bayes-optimal 
value of the MSE is quite large, and notably larger than the 
region in which the £1 minimization is able to give precise re- 
construction. Another point worth noting is that the location 
of the spinodal depends only very weakly to the value of the 
noise, for a large range of matrix and measurement noise (see 
inset of Fig. |2] and Fig. [3]): this shows that the robust from of 
the AMP algorithm is indeed robust to noise(s). 

6. CONCLUSIONS 

We have computed the Bayes-optimal value of the MSE for 
the reconstruction of sparse Gauss-Bernoulli signals in pres- 
ence of matrix uncertainty with the replica method, and con- 
sider a variant of the AMP algorithm robust to such uncer- 
tainty and to measurement noise. Finally, we have shown that 
AMP allows one to match the optimum MSE in a large region 
of parameters, and that the region is very weakly sensitive to 
measurement or matrix noises. Note that the present analysis 
applies to random i.i.d. measurement matrices; it is also pos- 
sible to use the seeded spatially coupled measurement matri- 
ces of |8, 13 , 14J and this would lead to an even larger region 
of optimality-matching performance for AMP. 
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